clear all
set more off

** RDD placebo FUA

global input "...\replication\data"
global output "...\replication\outputs"

use "$input\com_data.dta", clear

****************** Costmetics ******************
** Reomove Corsica
ren com_res com
drop if substr(com,2,1) == "A"
drop if substr(com,2,1) == "B"
destring com, replace

****************** Get dist border ******************
merge m:1 com using "$input\com_FUA.dta"
drop if _merge != 3
drop _merge
** Keep those within 50km of FUA border
drop if near_FUA_line < 0
drop dist_border
ren near_FUA_line dist_border
replace dist_border = -dist_border if near_FUA_poly != 0

****************** Get main vars******************
** Outcome
gen catch = ind if same == -1

egen sum_ind = sum(ind), by(com)
gen p_ind = ind/sum_ind

** Treatment
gen Treat = 0
replace Treat = 1 if same == -1

drop if dist_border == 0
destring dep_res dep_work, replace
*Due to proximity between FUA, drop mixed units
keep if dist_border < 6


****************************************************
** Com

rdbwselect p_ind dist_border, c(0) 
local opt_h=e(h_mserd)

rdrobust p_ind dist_border,c(0) h(`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_1

rdrobust p_ind dist_border,c(0) h(`opt_h')  covs(dep_res dep_work)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum p_ind if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_2

****************************************************
** Catchment
preserve 
keep if same == -1
rdbwselect catch dist_border, c(0) 
local opt_h=e(h_mserd)

rdrobust catch dist_border,c(0) h(`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum catch if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_3

rdrobust catch dist_border,c(0) h(`opt_h')  covs(dep_res dep_work)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum catch if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_4
restore

****************************************************
** Pop

collapse (mean) sum_ind dist_border, by(com dep_res)

rdbwselect sum_ind dist_border, c(0) 
local opt_h=e(h_mserd)

rdrobust sum_ind dist_border,c(0) h(`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum sum_ind if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_5

rdrobust sum_ind dist_border,c(0) h(`opt_h')  covs(dep_res )
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum sum_ind if abs(dist_border) < `opt_h' & dist_border < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_6

****************** Results in Table ******************
esttab ATE_1 ATE_2 ATE_3 ATE_4 ATE_5 ATE_6, tex ///
b(%9.3f) se(%9.3f) stats( h_l N N_bw m_outc, labels("Bandwidth" "Obs" "Obs. bw." "Mean dep. var. below" ) ///
fmt( %9.3gc %9.0gc %9.0gc %9.3gc)) nolines nogaps star(* 0.10 ** 0.05 *** 0.01) ///
keep(RD_Estimate) coef(RD\_Estimate  "ATE") 

****************** Results in Figure ******************
*Avoid bias due to departmental borders
drop if near_dep_b < 10

preserve
replace dist_border = -dist_border
rdbwselect p_ind dist_border, c(0) 
local opt_h=e(h_mserd)
rdplot p_ind dist_border if dist_border>-`opt_h' & dist_border<`opt_h', binselect(es) scale(2) p(2) ci(95) ///
graph_options(legend(off) scheme(s1color) title("Non-FUA ________________________________ FUA") ///
ytitle("Share of local residents working in FUA") ///
xtitle("Distance to the FUA border (km)")  legend(off))
graph export "$output\figure_placebo_FUA_noDep_vR2.pdf", replace
restore

preserve 
keep if same == -1
rdbwselect catch dist_border, c(0) 
local opt_h=e(h_mserd)
rdplot catch dist_border if dist_border>-`opt_h' & dist_border<`opt_h', binselect(es) scale(2)  p(2) ci(95) ///
graph_options(legend(off) scheme(s1color) title("Non-FUA ________________________________ FUA") ///
ytitle("Residential population") ///
xtitle("Distance to the FUA border (km)")  legend(off))
graph export "$output\figure_placebo_FUA_noDep_catch_vR2.pdf", replace
restore

rdbwselect sum_ind dist_border, c(0) 
local opt_h=e(h_mserd)
rdplot sum_ind dist_border if dist_border>-`opt_h' & dist_border<`opt_h', binselect(es) scale(2) p(2) ci(95) ///
graph_options(legend(off) scheme(s1color) title("Non-FUA ________________________________ FUA") ///
ytitle("Residential population") ///
xtitle("Distance to the FUA border (km)")  legend(off))
graph export "$output\figure_placebo_FUA_noDep_ind_vR2.pdf", replace